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ABSTRACT 

We present a new numerical code, PLUTO, for the solution of hypersonic flows in 1, 2 and 
3 spatial dimensions and different systems of coordinates. The code provides a multi-physics, 
multi-algorithm modular environment particularly oriented towards the treatment of astrophysi- 
cal flows in presence of discontinuities. Different hydrodynamic modules and algorithms may be 
independently selected to properly describe Newtonian, relativistic, MHD or relativistic MHD 
fluids. The modular structure exploits a general framework for integrating a system of conser- 
vation laws, built on modern Godunov-type shock-capturing schemes. Although a plethora of 
numerical methods has been successfully developed over the past two decades, the vast major- 
ity shares a common discretization recipe, involving three general steps: a piecewise polynomial 
reconstruction followed by the solution of Riemann problems at zone interfaces and a final evo- 
lution stage. We have checked and validated the code against several benchmarks available in 
literature. Test problems in 1, 2 and 3 dimensions are discussed. 



Subject headings: Methods: numerical - hydrodynamics - magnetohydrodynamics (MHD) 
tivity - shock waves 
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1. Introduction 

Theoretical models based on direct numerical simulations have unveiled a new way toward a better 
comprehension of the rich and complex phenomenology associated with astrophysical plasmas. 

Finite difference codes such as ZEUS rtStone fc Normanlll992allr} l or NIRVANA+ (|Zieglerlll998l ) inaugu- 
rated this nov el era and h ave been used by an increasingly large fraction of researchers nowadays. However, 
as reported in lFalkj ( 2002 ). the lack of upwinding techniques and conservation properties have progressively 
moved scientist's attention toward more accurate and robust methods. In this respect, the successful em- 
ployment of the so-called high resolution shock-capturing (HRSC) schemes have revealed a mighty tool to 
investigate fluid dynamics in non-linear regimes. Some of the motivations behind their growing popularity is 
the ability to model strongl y supersonic f lows while retain ing robustness and stability. The unfamiliar reader 
is refereed to the books of iTorol (jl997l ). iLeVeauel (|l998h and references therein for a more comprehensive 



overview. 
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Implementation of HRSC algorithms is based on a conse rvative for mulation of the fluid equations and 
proper upwinding requires an exact or approximate solution (Roc 1986) to the Riemann problem, i.e., the 
de cay of a d i sconti nuity separating two constant states. This approach dates back to the pioneering work 



Godunovl (|l959h . an d it has now becom e the lead line in developing high resolution codes examples of 



of 

which incl ude FLASH jFrvxell et alj |200cl for reactive hydrody namics), the special relativistic hyd ro code 
GEN ESIS |A!ov et alJll999h . the Versatile Advection Code (VAC. lT6th|[l996h or the new NIRVANA (|Ziegler 
l2004h . 

Most HRSC algorithms are based on the so-called reconstruct-solve-average (RSA) strategy. In this 
approach volume averages are first reconstructed using piecewise monotonic interpolants inside each com- 
putational cell. A Riemann problem is then solved at each interface with discontinuous left and right states 
and the solution is finally evolved in time. It turns out that this sequence of steps is quite general for many 
systems of conservation laws and, therefore, it provides a general framework under which we have devel- 
oped a multi-physics, multi-algorithm high resolution code, PLUTO. The code is particularly suitable for 
time-dependent, explicit computations of highly supersonic flows in presence of strong discontinuities and 
it can be employed under different regimes, i.e., classical, relativistic unmagnetized and magnetized flows. 
The code is structured in a modular way allowing a new module to be easily incorporated. This flexibility 
turns out to be quite important, since many aspects of computational fluid dynamics are still in rapid de- 
velopment. Besides, the advantage offered by a multi-physics, multi-solver code is also to supply the user 
with the most appropriate algorithms and, at the same time, provide inter-scheme comparison for a better 
verification of the simulation results. PLUTO is entirely written in the C programming language and can 
run on either single processor or parallel machines, the latter functionality being implemented through the 
Message Passing Interface (M PI) library. The code has already been succe ssfully employed i n the context of 



stellar and extragalac tic jets (jBodo et al 



Massaglia et al.l 120051 ). accretion disks (jBodo et al 



2003t iMignone et al.lbooi I2005I). radiative s hocks |Mignonell2o"ol 



2005 



Tevzadze et al.l Isubmittedl ). magneto-rotational 



instability, relativistic Kelvin-Hclmholtz instability and so forth. 

The paper is structured as follows: in <j2]we give a description of the code design; in f}3]we introduce the 
physics modules available in the code; in 2] we give a short overview on source terms and non-hyperbolicity 
and in S}5]the code is validated against several standard benchmarks. 



2. Code Design 

PLUTO is designed to integrate a general system of conservation laws which we write as 

^ = -V-T(U)+S(U). (1) 

Here U denotes a state vector of conservative quantities, T(U) is a rank-2 tensor, the rows of which are the 
fluxes of each component of U and S(U) defines the source terms. Additional source terms may implicitly 
arise when taking the divergence of T(U) in a curvilinear system of coordinates. An arbitrary number of 
advection equations may be added to the original conservation law (flj). 

Although the components of U are the primary variables being updated, fluxes are more conveniently 
computed using a different set of physical quantities which we take as the primitive vector V. This choice is 
supported, moreover, by the fact that interpolation on primitive variables enforces physical constraints such 
as pressure positivity and sub-luminal speeds in the case of relativistic flows. 

Numerical integration of the conservation law |T]) is achieved through shock-capturing schemes using 
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the finite volume (FV) formalism where volume averages evolve in time. Generally speaking, these methods 
are comprised of three steps: an interpolation routine followed by the solution of Riemann problems at zone 
edges and a final evolution stage. In PLUTO, this sequence of steps provides the necessary infra-structure 
of the code, see the schematic diagram in Fig. [5J 

At the higher level, the original system of equations is integrated by following this general sequence of 
steps, regardless any knowledge of the physics involved. The explicit form of U, V, T(U) and S(U), on the 
other hand, depends on the particular physical module selected. Thus, at the lower level, a physical module 
collects the set of algorithms required to compute the terms involved in the discretization of the right hand 
side of Eq. {!]). This set should provide one or more Riemann solver (s), mapper routines for the conversion 
between primitive and conservative variables, a flux routine giving the components of T(U) in each direction, 
a source term function (if any) and a routine to compute the maximum and minimum characteristic speeds 
of the Jacobian matrix. Of course, additional features may be easily added by exploiting the independent 
modularity. 

From the user's perspective, a particular configuration can be defined through a friendly interface entirely 
written in the Python scripting language. The interface allows the user to specify all problem-dependent 
attributes and algorithms, such as number of dimensions, geometry, physics module, reconstruction method, 
time stepping integration and so forth. 



2.1. Notations 



PLUTO employs logically rectangular grids in a generic system of orthogonal curvilinear coordinates 
(x , x 2 , x 3 ) (see Tabic [TJ. Let Ni, N 2 and N 3 be the number of points in the three directions. The lower 
and upper coordinate bounds of zone (i,j,k) _i) and j^+i ) x k+i)> res P ec tively, 



where 1 < i < Ni, 1 < j < Ns, 1 < k < N3. The zone widths are then simply given by Ax\ — x 1 L — x 1 



i+4 



Ax 2 — x 2 , 1 — x 2 1 , Ax?. = x, , 1 — x, 1 . The mesh spacing can be either uniform or (geometrical or 

3 3+2 3-2 h k +2 k ~2 y 
logarithmic) stretched. 

Parallelization is achieved through domain decomposition: the global domain is divided in sub-domains 
with an equal number of points and each of the sub-domains is assigned to a processor. By default the 
sub-domains are created as maximally cubic, however the user, at run time can specify a different strategy 
for the creation of the sub-domains, imposing that a given direction cannot be subdivided. In the code, 
parallelization is handled through the MPI library. 

Grid adaptivity technique s are also being provid ed. A one dimensional adaptive mesh refinement 
(AMR) module based on the IrBerger fc Colellal (| 19891) method is distributed with the code and exten- 



sion to multiple spatial dimensions is currently being developed using the CHOMBO library available at 



http://seesar.lbl.gov/ANAG/chombo/. However, since the main goal of this paper is to focus on the code 



modularity and its performance, AMR implementation will be described in future works. 

In order to avoid formulas with cluttered notations we omit integer-valued subscripts when referring to 
three-dimensional quantities. Thus ~Vi,j,k becomes simply V. We introduce the standard two-point difference 
one-dimensional operator 

£ d ( Y ) = -^d( A + F +- A - F -) + sd > ( 2 ) 

where d — 1, 2, 3 is a given direction, A± and AV d are, respectively, the cell's right (+) and left (— ) interface 
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Fig. 1. — Simplified flow diagram of the reconstruct-solve-average (RSA) strategy: first, volume averages 
U are more conveniently mapped into primitive quantities V. Left and right states V+.l and V-.r, are 
constructed inside each zone by suitable variable interpolation and/or extrapolation. A Riemann problem is 
then solved between V+,l and V+ i r to compute the numerical flux function F+ at cell interfaces and the 
solution is finally advanced in time. 



Table 1. Systems of coordinates (i.e. coordinates, volumes and areas) adopted in PLUTO and their 

meaning. 
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Note. — r is the cylindrical or spherical 
radius, < <^ < 27r is the azimuthal angle and 
< 8 < 7r is the polar angle. Here A n r = 
(r™ — r™)/n, A/i = cos 9— — cos 9+, s+ = 
sin 8+ , where + or — refer to the right or left 
zone interface, respectively. 
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areas and cell volume in that direction (see Table [I]). Here ± = (i ± ^, j, k), ± 5, fc), (i, j, fc ± 5) when 
(i = 1,2,3, respectively. Scalar components of the vectors appearing in Eq. ^ will be denoted with a 
pair of square brackets; e.g. £^ m0 ] * s * ne component of £ contributing to the equation. Furthermore, 
since several of the building block algorithms are one-dimensional, we will also omit the superscript d when 
unnecessary. 

The numerical flux functions F± in Eq. follow the solution of 1-D Riemann problems at cell 
interfaces. Since more than one Riemann solver may be available in each physics module, we set, without 
loss of generality, 

F+=^(V + , L ,V +iR ) , (3) 

where 1Z is the Riemann solver being used and V+x, V+,r are suitable left and right states at the zone 
edges perpendicular to the x direction. 



2.2. Reconstruction 



Interpolation routines are designed to reconstruct a piecewise polynomial approximations "P(x) to V 
inside each cell starting from its cell averages: 



V± lS = T(P,V) , 



(4) 



where S = L (S = R) at x = x + (x = x_) . Here I is an interpolation routine designated to provide left and 
right edge interpolated values inside each cell, that is, V +! l = lim x _> x+ 'P(x) and V_r = lim x _ +x _ 'P(x), 
where the "L" and "R" subscripts refer to the sides of the interface. 

Reconstruction methods have to satisfy monotonicity constraints in order t o avo i d spurious oscilla tions 
in proximity of discontinuities and steep gradients, see the books from iTorol |l997t ). iLeVequel |l998h and 
references therein. For second-order linear interpolants, for instance, one has 



V ±iS = V± 



AV 



(5) 



with the slopes AV computed following a limiting procedure applied to primitive or characteristic variables. 
In the latter case one has 



AV 



w k r k , Aw k = ]im(Au>k i+1 Aw k ,-) 



(6) 



where r k and If. are, respectively, the right and left eigenvectors of the primitive form of the equations, 
Awk,± = ±lfc • (Vj±i — Vj) are forward (+) and backward (— ) derivatives and k labels the fc-th charac- 
teristic field. Different slope limiters (lim) are characterized by distinct steepening properties and can be 
independently assigned to each primitive variable or characteristic field. 

At the time of this writing, PLUT O allows to choose between eit her flat (1 st order in space), lin- 
ear (2 nd order), thir d-order convex ENO (jDel Zanna fc Bucciantinil 120021 ) , parabolic reconstr uctions (as in 
Mignone et al.ll2005T ). or the 5th order finite difference WENO scheme of I Jiang fc Shul (|l996h (only for the 
hydrodynamics equations) . 
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2.3. Riemann Solver 



As already anticipated, computation of the numerical flux function F + at a zone edge (#+) requires the 
solution \J{x, t), for t > to the initial value problem 



U(a,i ) 




if 
if 



X < x + , 
X > x + , 



(7) 



complemented with the one-dimensional evolutionary equation for U. Since this section deals with interface 
quantities, we shall omit, in this section only, the + subscripts. 

The exact solution to the Riemann problem involves the decay of a set of non-linear waves and can 
be a rather cumbersome task to achieve. With the exception of few simple cases, existing Riemann solvers 
routinely involved in upwind schemes are based on different levels of approximation. 

The Lax-Friedrichs Rusanov flux is robust, but also the most diffusive solver and is available for all 
modules. It computes the fluxes according to: 

„ 1 



f L + fa- |A max | (U R - U L ) 



where f = e d ■ T(U) is the projection of the tensor flux on the e d 
Kronecker-Delta symbol and |A max | is the largest local signal velocity. 



(8) 

{Sid, fod, Ssd) unit vector, Sij is 



A somewhat different approximation comes from the HLL-family solvers: 

f L if Ai > , 

f fc if A fc A fe+1 <0 (fc = l,...,JV-l), 
f R if A w < , 



(9) 



where the solution to the Riemann problem is approximated by TV waves with Xk+i > A& and k = 1, . . . , N— 1 
and separated by N + 1 states. The ffc are computed from a suitable "parameterization" of the Rankinc- 
Hugoniot jump conditions across each wave: 



Afe (Ufe + i — Ufe) — ffc + i — ik 



(10) 



with Ujv+i = Ur. The HLL and HLLC solvers can be consistently d erived following this recipe w ith N = 2 
and N = 3, respecti vely, see Toro et al. ( 1994 ). Batten et al. (1997) for the hydro equations. O ( 20051 ) for 
the MHD equations, iMignone fc Bodol (|200.4 l200d ) fo r the relativistic equatio ns. Similarly, we have also 
implemented the multi-state (iV = 5) HLLD solver of ( Mivoshi fc KusancT 2005) in the MHD module. The 
HLL approach has some attractive features in that it guarantees positivity of pressure and, in the case of 
relativistic flows, it preserves the condition |v| < 1. 

Linearized Riemann solvers are more accurate since the averaging process inherent to Eq. ([8]) or (|9|) is 
removed and all jumps are computed by linearization around some average state. The Roe solver computes 
the numerical fluxes according to: 



f L + f R - J2 l A *l L fc • ( u l - u r) R fc 



(ii) 



where the rows (columns) of L (R) are the left (right) eigenvectors of the Jacobian df{\J)/d\J. 
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The Advection Upstream Splitting Method (AUSM, originally proposed Liou fc Steffen 1993h provides 
an alternative and efficient upwinding strategy by splitting the flux into convective and pressure terms, 
respectively associated with linear and nonlinear fields: 



F = &v n 



P, 



(12) 



where v n is a suitably defined convective velocity and the p flux is g overned by the a coust ic speed. The 
origin al AUSM scheme has been substantially improved in the work of lLioul j 19961 120061) and lWada k. Lioul 
|l997h . 



The most accurate approach consists in directly solving the whole set of Rankine-Hugoniot jump con- 
ditions to find V*, from which the flux can be computed, 



F = f(V*). 



(13) 



However, this approach is computationally expensive since it generally involves the solution of highly non- 
linear equations. 

Different physics modules come with different sets of Riemann solvers, and additional methods of solution 
may be easily accommodated. We warn that more accurate Riemann solvers (such linearized or exact 
schemes) may introduce insufficient numerical dissipation for certain flow configurations. Sporadically, this 
could lead to a number of numerical pathologies such as the c arbun cle phenomena, odd-even coupling or 
lack of positivity- preserving properties. See the work bv lQuirkl (|19941 ) for a comprehensive review. 



2.4. Temporal Evolution 

Time marching algorithms provide a general (i.e. physics-independent) framework for the discretization 
of the left hand side of Eq (fT|). For example, in the simplest case of forward Euler discretization one has 

U" +1 - U™ 

— sr-=* B . (") 

where C n is the right hand side operator ([2]) or a sum of them for dimensionally s plit or unsplit schemes 



respe ctively. The time step At is limited by the Courant-Friedrichs-Lewy (CFL. ICourant. Friedrichs fc Lewv 
19281 ) condition: 

/ Al d ■ \ 

Af = C a min -pL , (15) 

d VI A maxl/ 

with A^ in and A^ ax being, respectively, the smallest cell length and largest signal velocity in the d direction. 
C a identifies the Courant number, restricted by the conditions given in the last column of Table [2] 

PLUTO provides a number of time-marching schemes for the explicit numerical integration of the 
conservation law (fTJ). We discriminate between 1) fully discrete, zone-edge extrapolated and 2) semi-discrete 



methods. Evolution in more than one dimension may be achieved by either operator splitting (|Stranelll968f ) 
or fully multidimensional integration. 
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2.4-1- Zone-Edge Extrapolated Methods 



Zone-edge extrapolated methods achieve second order temporal accuracy by midpoint in time quadra- 



ture: 



U n+1 = U" + At£ n+ i . (16) 
and thus are based on a single-step. For directional splitting methods one has C n+ ^ = C d ^V n+ 5^ , whereas 
= Ed £d (v"+5) in the case of a dimensionally unsplit method. The input states for the Riemann 



(17) 



solver are computed by Taylor expansion of the primitive variables at half time step, 

At 



V" + ^ = V" ± - I I — ^A' 
1 Ax 



AV r ' 



where A" is the matrix of the quasi-linear form of the equations and I is the identity matrix. An alternative 
form which does not require the primitive form of the equations can be written in terms of conservative 
variables. Both variants yield the well-known MUSCL-Hancock scheme (MH henceforth, Ivan Leerl 11974 ; 



Torolll997f ). When eigenvalues and eigenvectors are available, upwind limiting may be used to select only 
those characteristics which contribute to the effective left and right states . Thi s approach is employed, 
for instance, in the original PPM advection scheme of Colella 8z Woodardl (1984) and we will refer to as 
characteristic tracing (ChTr henceforth). These methods require boundary conditions to be assigned at the 
beginning of the time step only. 

The advantage offered by dimensionally split single-step algorithms is the lower computational cost (only 
one Riemann solver per cell per direction is required, see Tab[2|) and storage requirement. An alternative, 
di mensionally un s plit versio n of th ese schemes, however, adopts the corner transport upwind (CTU) method 
of lColellal (|l990h . ISaltzmanl (|1994l ) and is more expensive. In this case, an extra correction term is needed 
in Eq. (| 1 7[) : in two dimensions, for example, the input states for the Riemann problem are modified to 



fr" + 3 
u ±,s 



U 



±,S 



At f t,n + i 



(18) 



with £*> n +3 being the right-hand side operator corresponding to the transverse direction. As in the dimen- 
sionally split case, the CFL stability condition is not affected by the dimensionality of the problem (i.e. 
C a < 1, see also Tab [2]). Contrary to its dimensionally split version, on the other hand, this method is 
computationally more expensive since 4 instead of 2 (in 2D) and 12 instead of 3 (in 3D) Riemann problems 
must be solved at each interface, see Table [2] 



2-4-2. Semi-Discrete Methods 



Semi-discrete methods are based on the classical method of lines, where the spatial discretization is 
considered separately from the temporal evolution which is left continuous in time. In this framework Eq. 
(fTJ) is discretized as a regular OD E. PLUTO impl e ments the 2 nd and 3 rd order Total Variation Diminishing 
(TVD) Runge-Kutta schemes of Gottlieb fc Shu ( 19961 ). The second-order method (RK2) advances the 
system of conservation law JT]) as 

U* = U" + AtC n , (19) 
1 



U n+1 = - 



U" + U* + AtC 



(20) 
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with C n = C d (V n ) or C n = ^2d,£ d (V n ) in the case of a dimensionally split or unsplit method, respectively. 
The third-order Runge-Kutta method (RK3) may also be used, at the cost of an additional step: 



U* 



U 



n+l 



U* = IT 



1 

= 4 
1 
3 



3U" + U* + AtC* 



U n + 2U* 



2At£* 



(21) 
(22) 



(23) 



For this class of methods, input states for the Riemann solver are given by the output of the interpolation 
routine, see N2.2I Besides, boundary conditions must be assigned before each step. A total of two and three 
Riemann problems per cell per direction must be solved by the RK2 and RK3 marching scheme, respectively. 
Furthermore, fully unsplit Runge-Kutta integrators require a stronger time step limitation, see Table [2j 



3. Physics Modules 

PLUTO is distributed with four independent physics modules for the explicit numerical integration of the 
fluid equations under different regimes and conditions. The hydrodynamics (HD), magnetohydrodynamics 
(MHD), relativistic (RHD), and relativistic MHD (RMHD) modules solve, respectively, the Euler equations 
of gas dynamics, the ideal/resistive MHD equations the energy-momentum conservation laws of a special 
relativistic perfect gas, and the equations for a (special) relativistic magnetized ideal plasma. 

In what follows, p, p, and E will denote, respectively, the proper density, thermal pressure and total 
energy density. Vector fields such as m = (mi, m 2 , m 3 ) T , v = (vi, 1)%, v 3 ) T , B = {B\, B 2 , B 3 ) T , define the 
momentum density, velocity and magnetic field. Finally, T will be used to define the ratio of specific heats 
for the ideal equation of state. 



3.1. The Hydrodynamics (HD) Module 

This module implements the equations of classical fluid dynamics with an ideal equation of state. The 
conservative variables U and the flux tensor are: 



(24) 



where m = pw is the momentum density and I is the unit, 3x3 tensor. The total energy density E is related 
to the gas pressure p by the ideal gas closure: 





f p \ 




( " V \ 


u = 


m 


, T(U) = 


mv + pi 




I E J 




\ (E+p)v J 



E 



P 



(25) 



r - 1 2 P 

The set of primitive variables V = (p,v,p) T is given by density, velocity v and thermal pressure p. 

This module comes with a set of several Riemann solvers, including the no nlinear Riemann solver base d 
on the two-shock approximations ( Colella fc WoodardI 1984 : Frvxell et al. 2000( l. the Roe solver ( Toro 1997 ). 
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the AUSM+ scheme 
Lax-Friedrichs solver (|Rusanovlll96l[ ) 



jLiorJll996h, th e HLL (jEinfeldt et alJll99l|) . HLLC (jToro et alJll994h solvers and the 



The HD module co ntains an imp lementation of the fast Eulerian transport algorithm for differentially 
rotating disk (FARGO. iMassetl 120001 ) on polar grids. The FARGO scheme allows much larger time steps 
than the standard integration where the Courant condition is traditionally limited by the fast orbital motion 
at the inner boundary. 



3.2. The Magnetohydrodynamics (MHD) Module 



The MHD module deals with the equations of classical ideal or resistive magnetohydrodynamics (MHD). 
In the ideal case, U and T may be written as 

. T 



U = 



/ 9 \ 
m 

B 

V E J 



T(U) = 



pv 

mv BB +p t I 
vB - Bv 

V (B + ft )v-(vB)B / 



(26) 



with m = pv andpt = P+|B| 2 /2 being the total (thermal + magnetic) pressure, respectively. The additional 
constraint VB = complements the magnetic field evolution (see fc|3.2.ip . Resistivity is introduced by adding 
appropriate parabolic terms to the induction and energy equations, see Sj4] 



Available equations of state implemented arc the ideal gas law, 

12 



E = 



P 



r - 1 



m • 



IBI 



(27) 



and the isothermal equation of state p = c 2 s p where c s is the (constant) isothermal speed of sound. 



The set of primitive variables is the same one used for the HD module, with the addi tion of magnetic 
fields. The user c an choose among the followi ng ava ilable Ri emann solvers: the Roe so lver of lCargo fc Gallice 
(jl997h . the HLL |janhunej2000[ ). HLLC |l1 |2005| ). HLLD (jMivoshi fc Kusanolliooj) and the Lax-Friedrichs 
solvers. 



3.2.1. Solenoidal Constraint 



The solution to the MHD equations must fulfill the solenoidal constraint, V ■ B = 0, at all times. Unfor- 
tunately it is well known that numerical scheme do not naturally preserve this condition unless special dis- 
cretization techniques a re used. Among the variety of monopole control stra tegies proposed in literature (for 



a review see 



TotbfeOQOl ). we have implemented 1) the ei ght wave formulat ion dPowelJl994HPowell et alfl 999) 



and 2 ) the constrained transport (CT henceforth) of iBalsara fc Spicerl (J1999J), and lLondrillo fc Del Zanna 



(200 4). The CT fram e work has been incorporated into the u nsplit CTU integrato r following the recent work 
by Gardiner fc Stone ( 2005 ). A similar approach is used by Tevssier et al. ( 20061) . 



In the first strategy, the magnetic field has a cell centered representation and an additional source term 
is added to the MHD equation. T he discretization of the source term is different depending on the Riemann 
solver (following Janhunenl 2000 ). a feature which we found to greatly improve robustness. 



In the CT formulation, on the other hand, the induction equation is integrated directly using Stoke's 
theorem and the magnetic field has a staggered collocation. 



3.3. The Relativistic Hydrodynamics Module (RHD) 



This module deals with the motion of an ideal fluid in specia l relativity. The equation s are given by 
particle number conservation and energy-momentum conservation ( Landau fc Lifshitall959 ). Conservative 
variables and tensor flux are: 





( D \ 




/ Dv \ 


u = 


m 


, T(U) = 


mv + pi 




V E ) 




V m / 



(28) 



where D = jp is the laboratory density and 7 = (1 - 
are normalized to the the speed of light. 



|2\-i 



2 is the Lorentz factor. For convenience, velocities 



The relation between conserved and primitive variables V = (p, v,p) T is more involved than its classical 
counterpart. Specifically we have 



D = 7/9 , m = phj 2 v , E — ph^ 2 



P, 



(29) 



where the specific enthalpy h = h(p, p) depends on the equation of state. The inverse map requires the 
solution of a nonlinear equation. iRvu et al. address this topic for different choices of h(p, p). 



Two equations of state are available for this module: the constant-I" ideal gas law 

r 



f 



f 



-6. 



with O = p/p and the equation of state (TM henceforth) already introduced in lMignone et al. I (|2005h : 



h = -e 
2 



9 

-e 2 

4 



1, 



(30) 



(31) 



which satisfies Taub's inequality (|Taublll948) and approximates (within < 4%) the single-specie relativistic 
perfect gas (for a thorough discussion, see Mignone et al. 2005 : Rvu et al. 20061 and references therein). 



The two shock nonlinear Riemann solver described in lMignone et al.l (|2005l ) has been i ncorporated into 
the se t of available Riemann solvers, together with the recently developed HLLC algorithm bv lMignone fc Bodol 
(|2005l ). The HLL and Lax-Friedrichs Riemann solvers have also been coded. 



3.4. The Relativistic Magnetohydrodynamics Module (RMHD) 



The motion of an ideal relativistic magnetized fluid is described by conse r vation o f mas s, energy- 
momentum, and by Maxwell's equations, see, for example, Unile fc Pennisil |l987h : Unile I (| 19891) . PLUTO 
implements the equation of special relativistic magnetohydrodynamics where the vector of conservative vari- 
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ables and respective fluxes can be cast as 

/ D \ 



U 



m 
B 



T(U) 



WtJ vv 

vB 



Dw \ 
- bb + Ip t 
Bv 

/ 



(32) 



where w t = ph + 6 2 n , &^ = |B| 2 /7 2 + (v • B) 2 , b = B/7 + 7(v • B)v and p t = p + 6^/2 is the total (thermal 
+ magnetic) pressure. 

The components of U are related to the primitive variables V = (ft, v,p, B) through 

D = pj, 

m = (^ 7 2 + |B| 2 )v-(vB)B, 

|B| 2 , |v| 2 |B| 2 -(vB) 2 



E 



phj — p 



(33) 
(34) 

(35) 



The inverse map is highly nonlinear and different methods of solution have been proposed when a constant 
r law is adopted (fo r a rec ent review, see iNoble et all l2006h . In PLUTO, we follow the approach described 



Mignone fc Bodd ([20061 ) which was recently generalized t o include the TM equat i on of state described in 



33.31 The interested reader should also refer to the work by iMignone fc Mc Kinnevl (|2007t ) . 

The RMHD module may be us ed with the Lax-Friedrichs, HLL or the recently developed HLLC Riemann 
solver, see Mignone fc Bodo ( 20061 ). Since the induction equation takes the same form as in the classical case, 
the divergence of the magnetic field can be kept under control by any of the methods already introduced in 
the MHD module, see 33.2.11 



4. Source Terms and non-Hyperbolicity 

In addition to the homogeneous terms previously described, a number of additional features may be 
added in the code. 

Local source terms are functions of the variables themselves but not of their derivatives and are included 
either during the advection step or via operator splitting. Examples include the centripetal and Coriolis terms 
implicitly arising when working in polar or spherical coordinates, external forces (e.g. gravity) or optically 
thin radiative losses. 

Non-ideal effects such as viscosity, resistivity and thermal conduction, on the other hand, introduce 
parabolic corrections to the equations and involve the solution of diffusion equations. 



4.1. Geometrical Source Terms 



The divergence of a rank-2 tensor in curvilinear coordinates breaks down the homogeneous properties of 
Eq. ([TJ) . This loss of conservation comes from the additional source terms inevitably in troduced by those uni t 
vectors that do not have fixed orientation in space (see, for instance, the Appendix in Mignone et al. 2005h . 
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In the case of symmetric or anti-symmetric tensors, however, some of the source terms can be eliminated 
by properly re-arranging the derivatives. This ensures conservation of the angular momenta rather than the 
linear ones. 

For all physics modules previously introduced, for example, the components of the flux tensor in the 
momentum equation form a 3 x 3 rank 2 symmetric tensor My +p6ij = Mji + pSji, where p is the isotropic 
pressure term. However, since V • (Ip) = Vf>, the differencing terms involving pressure are separately treated 
as gradient components and never contribute to the source terms. 

The symmetric character of M leads to further simplifications in polar coordinates, since the (^-component 
of the divergence of M may be written as 

(V .M V = 4^) + I^ + ^. (36) 
r z or r oq> oz 

This leads to the so-called "angular momentum-conserving form" where the radial contribution to is 
computed by the new operator 



1 



rAV 



riFl, m ,-riFL, m A , (37) 



thus leaving the radial component of momentum as the only non-homogeneous equation, with Sr mr ] = M^/r. 
Similar considerations lead to the angular momentum conserving form in spherical geometry. 

Geometrical source terms are integrated explicitly during the advection step by adding their contribution 
to the right hand side operator ©. In the semi-discrete approach, this is straightforward. For edge- 
extrapolated methods, this requires 1) augmenting (|17p with the cell-center source term for a half step 
and then 2) averaging the resulting time centered left and right edge values to compute the S for the final 
conservative step. 



4.2. Optically thin cooling 



For many astrophysical applications, radiative processes may become important during the evolution. 
This is particularly true when the cooling timescale becomes comparable to, or smaller than the typical 
timescale for the dynamical evolution. 

Optically thin radiative losses are treated as local source terms depending on the temperature, density 
and an arbitrary number of ions of different elements (e.g. H, He, C, N, O, and so on). The ionization 
fractions are advected with the fluid and are subject to collisional ionization, recombination and charge 
exchange processes. 



Curr ently, we have imp lemented a cooling module (jTesileanu et aLllin preparation! ) similar to the one 
coded by iRaga et al.l (|1997l ). extending the temperature range of applicability (2000 < T < 2 • 10 5 K) by 
using an increased number of ion species. This module employs 28 ions and is valid for densities up to 10 5 
cm~ 3 . 

Operator splitting is employed to evolve the chemical network in time. A 2 nd order fully implicit scheme 
is adopted in regions of rapid variations, whereas a second order explicit integrator is used otherwise. A 
comprehensive description is outside the goal of this work and will be available in future works. 
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4.3. Treatment of Parabolic Terms 



Parabolic terms introduce second-order spatial derivatives and their treatment requires the solution of 
diffusion equations. Typical examples include electric resistivity, 



<9B 

— + V x (-v x B + 77J) = 0, 

— +V-[(.E+p t )v-(vB)B + » 7 JxB] = 0, 
where J = V x B is the current and 77 is the resistivity; or thermal conduction, 



dE 
~dt 



V • (F 



adv 



kVT) 



0. 



(38) 
(39) 

(40) 



where F adv is the energy advection flux and k is the thermal conductivity coefficient. These terms may 
be included into the original conservation law with the further time step limitation At = min(At ad , Ai par ), 
where Ai ad is the advective time step (Eq. [T5|) and 

x 2 " 

(41) 



Ai par < 0.5 min 

d=l,2,3 



(Ax d ) 



max(a) 



is the diffusion step. Here a characterizes the physical process, e.g., r\ or k. 

For advection-dominated problems (i.e. Ai ad > Ai par ), Eq. pij) does not pose severe restrictions and 
diffusion terms can be treated explicitly. We achieve this by adding centered in space, second-order finite 
difference approximations to the right hand side operator £ d , see Eq. ©. 

However, for diffusion-dominated problems and/or increasingly high resolution simulations, At par will 
eventually drop below the typical advection scale. In such situa tions the explicit integration can be consid- 
erably accelerated using the super time stepping (STS) method ( Alexiades et al. 19961 O'Sullivan fc Downed 
2006). In STS, diffusion terms are included via operator splitting, and the solution vector is evolved over 
a super time step AT consisting of N smaller sub-ste ps, the stability of w hich is closely related to the 
properties of Chebyshev polynomials. It can be proved ( Alexiades et al. 19961) that 



AT -► N 2 At pa 



(42) 



which makes STS almost N times faster than the standard explicit scheme. Thus, if AT is taken to be 
the advective time step, STS will require (approximately) ^/Ai ad /At par iterations rather than At ad /Ai par , 
typical of a normal sub-cycling explicit time stepping. This approach offers dramatic ease of implementation 
over implicit schemes. 

Resistivity has been extensively tested through the direct comparison with analytical solutions of linear 
diffusion problems of magnetic field. Our results confirm the efficiency and accuracy of the STS approach 
already highlighted by previous investigators. This has encouraged further experiments towards nonlinear 
problems (e.g. Spitzer-like conductivity), where preliminary results suggest that the STS methodology may 
be succesfully applied. This issue will be presented in future works. 



Code Verification 



PLUTO has been validated against several be nchmarks typically ad opte d for other numerical schemes . 
Several algorithms for relativistic flows presented in iMignone et all (|2005l) and lMignone & Bodol |2005ll2006l ) 
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Table 2. Number of Riemann problems per cell per time step required by different time marching schemes 

in Nd = 1,2,3 dimensions. 



Time Marching 


1-D 


2-D 


3-D 


£*max 


MH (S) 


1 


2 


3 




i 


MH (U) 




4 


12 




i 


ChTr (S) 


1 


2 


3 




i 


ChTr (U) 




4 


12 




i 


RK2 (S) 


2 


4 


6 




it 


RK2 (U) 


3 


6 


9 


V- 




RK3 (S) 


2 


4 


6 




i 


RK3 (U) 


3 


6 


9 


V- 





Note. — (S) or (U) stand for dimensionally split 
or unsplit algorithm, MH and ChTr are the MUSCL- 
Hancock and Characteristic tracing schemes, RK2 
and RK3 are the 2 nd and 3 rd order Runge-Kutta 
schemes. The rightmost column gives the maximum 
allowed Courant number. 

+ High order interpolants (PPM, CEN03) may re- 
quire a lower limit. 



Table 3. Numerical schemes adopted in the different test problems described in the text. 



Test 


Module 


Case 


Geometry 


Dim 


Interpolation 


Solver 


Time Stepping 


C a 


Double Mach ReflectionfflBTt 


HD 


(a) 


Cartesian 


2 


Parabolic 


HLLC 


RK3 (S) 


0.8 






0) 






CEN03 


HLLC 


RK3 (S) 


0.8 






(c) 






WEN05 


Roc F-S 


RK3 (S) 


0.8 






(d) 






Parabolic 


Two-Shock 


ChTr (S) 


0.8 






(e) 






Linear 


Roc 


MH (U) 


0.8 






(f) 






Linear 


HLL 


MH (U) 


0.8 


Undcr-cxpandcd Jet ( ^5*2l 


HD 


(a) 


Cylindrical 


2 


Linear 


AUSM+ 


MH (U) 


0.4 






(b) 






CEN03 


HLL 


RK3 (S) 


0.4 


Rotated Shock Tube (i|5.3[ 


MHD 


(a) 


Cartesian 


2.5 


Linear 


Roe 


RK2 (U) 


0.4 






0) 






Linear 


HLLD 


RK2 (U) 


0.4 






(c) 






Linear 


HLL 


RK2 (U) 


0.4 


Fast Rotor (H5.4I 


MHD 


(a) 


Cartesian 


2 


Parabolic 


HLLC 


ChTr (U) 


0.6 






(b) 


Polar 




Linear 


Roc 


RK3 (U) 


0.4 


Magnetized Torus (H5.51 


MHD 


(a) 


Spherical 


2.5 


Linear 


HLLD 


R.K2 (U) 


0.4 






0) 


Spherical 




Parabolic 


HLL 


RK3 (U) 


0.4 






(c) 


Cylindrical 




Linear 


HLLD 


RK2 (U) 


0.4 






(d) 


Cylindrical 




Parabolic 


HLL 


RK3 (U) 


0.4 


Spherical Shock Tube (i|5.5t 


RHD 


(a) 


Spherical 


1 


Parabolic 


Two-Shock 


ChTr(S) 


0.8 






(b) 


Spherical 


1 


Linear 


HLL 


MH (S) 


0.8 






(c) 


Cartesian 


3 


Linear 


HLLC 


ChTr(S) 


0.8 


Magnetized Blast Wave (H5.7I 


RMHD 


(a) 


Cartesian 


3 


Linear 


HLL 


R.K2 (U) 


0.2 






(b) 


Cylindrical 


2 


Linear 


HLL 


MH (U) 


0.4 



Note. — Here ChTr stands for the characteristic tracing, RK2 and RK3 for the 2 nd and 3 rd Runge-Kutta time stepping, 
whereas MH is the MUSCL-Hancock time marching scheme. Split or unsplit schemes are denoted with (S) or (U). The rightmost 
column gives the CFL number C a . 
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are now part of the code and the related verification problems will not be repeated here. In what follows, 
we propose additional tests aimed to check the code performances on problems of different nature, geometry 
and dimension. Computational details for each test are given in Table [3] 



5.1. Double Mach Reflection of a Strong Shock 



The initial condition for this test problem (j Woodward fc Coleilalll984f) consists in a planar shock front 



making an angle a — 7r/3 with a reflecting wall, taken to be the x axis: 



(l.4,0,0,l) for x>x s (0) 

(p,v x ,v y: p) = < / \ (43) 

v ' I 18,8.25,-8.25,116.5) otherwise, 



where x s (t) — (lOt/ s'ma + 1/6 + yj tana) is the shock coordinate. The ideal equation of state with T = 1.4 
is used throughout the calculation. The computational domain is the rectangle [0,4] x [0, 1] covered with 
a uniform grid with 1/Ax — 1/Ay = 480. Outflow boundary conditions apply at x — 4 and reflective 
conditions are imposed at the bottom boundary y = for x > 1/6. Fluid variables are kept constant to their 
initial values at y — for x < 1/6 and at the leftmost boundary x — 0. At the top boundary y = 1 the exact 
motion of the shock is prescribed. We carry out integration until t = 0.4 using the six different combinations 
of algorithms described in Table [3l Panels in the left column of Fig. [2] show, from top to bottom, the results 
obtained with high order (> 2 ) interpolants : piece wise parabolic reconstruction (a), 3 rd order central ENO 
(b ) and the WENO scheme of Jiang fc Shu ( 19961 ) (c). Panels on the right show the original PPM scheme 



of IColella fc Woodardl (|1984h (case d, top), the 2 nd order unsplit Muscl-Hancock with a Roe solver (case e, 



middle) and with the HLL solver (case f, bottom). 

After the reflection, a complicated flow structure develops with two curved reflected shocks propagating 
at directions almost orthogonal to each other and a tangential discontinuity separating them. At the wall, a 
pressure gradient sets up a denser fluid jet propagating along the wall. Kelvin- Helmholtz instability patterns 
may be identified with the "rolls" developing at the slip line. This feature is visible only for some of the 
selected schemes. Indeed, the use of high-order interpolants, such as PPM (a,d) or WENO(c), and the 
ability of the Riemann solvers to resolve contact and shear waves result in smaller numerical viscosity when 
compared to a second-order slope-limited reconstruction and/or a more approximate Riemann solver. In 
this respect, case a) shows the smallest amount of dissipation, whereas the HLL solver (f, bottom right) has 
the largest numerical diffusion. 

In terms of CPU time, we found that case a) to f) followed the ratios 1.82 : 1.55 : 8.52 : 0.94 : 1.24 : 1, 
that is, with the 5 th order WENO schemes being by far the most expensive (more than a factor of 8 compared 
to the simple HLL) and the original PPM (split) scheme being the fastest. In doing the comparison, however, 
one has to bear in mind that the amount of computing time depends, in the first place, by the number of 
Riemann problems solved in each zone at each time step, as shown in Table [2] This number is 6 for cases a), 
b) and c), but 4 for case e) and f) (which are dimensionally unsplit) and only 2 for the original dimensionally 
split PPM scheme, case d). Thus one can safely conclude, comparing case e) and f), that the Roe solver is by 
a factor ~ 1/4 slower than HLL or, by comparing a), b) and c), that the CENO and PPM are considerably 
faster than WENO, although the combinations of algorithms in a) better resolves small scale structures. 

Finally, we note that computations carried with more accurate schemes (with the exception of the 5 th 
order WENO) exhibit a slight tendency to form a spurious "kinked" Mach stem on the x axis. This feature 
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Fig. 2. — Density maps for the double Mach reflection test at t = 0.2. Each panel shows results obtained with the different 
combinations of schemes listed in table [3] The mesh size (1/Ax = 1/Ay = 480) and Courant number (C a = 0.8) are the same 
for all cases. For the sake of clarity only the region [0, 3] X [0, 1] is shown. 
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is engendered by a numerical flaw (|Quirklll994[ ) caused by insufficient numerical dissipation and it becomes 
particularly dramatic when the resolution is further increased. 



5.2. Under-expanded Jet 

Code validation, in addition to the typical tests described above, can be carried out against laboratory 
experiments of under-expanded free jets as well. This kind of experiments consists of injecting, through a 
converging nozzle, a gas at stagnation pressure po into a chamber kept at lower pressure p c . The shock 
structure that forms is typically axially symmetric and consists of a quasi-statio nary Mach-disk, located at a 
distance z m from the nozzle, and of a barrel and reflected shocks ( Youne 1975 ) . After series of experiments 



with different gases, an em pirical expression was found ( Adamson fc Nicholls 1959t Bier fc Schmidt 



Ashkenas fc Sherm an 1966) that relates the Mach-disk location z m with the pressure ratio (jYoune 



197 



1961 

ST 




(44) 



with r* the effective sonic nozzle radius. Furthermore, iKnuthl (|1964I ) obtained for the Mach number M a 
on the jet axis the following empirical expression: 

r+i 



M axls «(2.2)— r(r-i 



i 



r - i 



Gr*) 



r-i 



(45) 



We have carried out 2D numerical simulations in cylindrical coordinates on a uniform grid r £ [0,40], 
z 6 [0, 80] with 20 zones per beam radius. The grid has been further extended in both directions (up to 
r = 80 and z — 160) by adding a second patch of geometrically stretched zones (100 and 200 zones in r 
and z, respectively) in order to avoid spurious reflections at the boundaries. Free outflow is set at the outer 
boundaries and reflective conditions are imposed on the axis r = and for r > 1 at z = 0. Considering 
the actual nozzle as the injection zone at z = 0, r < 1, we hav e obtained the values of the pressure p* and 
density p* by employing the isentropic laws for a perfect gas ( Shames! 1983 ) for a converging nozzle with 
stagnation pressure po and density po: 



Po 



r + i 



Po 



r + i 



(46) 



We choose po jp c = 2 ■ 10 3 and po/p c — 2 • 10 4 (typical of an argon jet in helium at ambient temperature) and 
accordingly choose a sonic injection velocity at the the nozzle (M* = 1). Densities, velocities and lengths 
are normalized to the ambient p, sound speed and beam radius, respectively. 

Fig. [3] shows the density maps for the two selected numerical schemes (see Table [3]) at t = 240. The 
results obtained with the AUSM+ Riemann solver (scheme a) disclose a somewhat greater amount of small 
scale structure than the fFLL integration (scheme b). This is not surprising because of the former's ability to 
properly capture contact and shear waves. On the contrary, the latter greatly simplifies the wave pattern of 
the Riemann fan at the cost of introducing extra numerical dissipation. This is compensated, to some level, 
by the choice of the 3 rd order interpolant (CEN03) which, nevertheless, is better combined with an equally 
accurate time marching scheme (RK3). The overall balance, therefore, favors the second order scheme since 
only 4 Riemann problems must be solved instead of the 6 required by RK3, see Table [5J This conclusion 
is supported not only by the decreased numerical viscosity inherent to scheme a), but also by the reduced 
computational cost, which sees scheme b) being ~ 1.6 times slower than scheme a). 
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Fig. 3. — Density logarithms for the under-expanded jet at t = 240, for the linear Muscl-Hancock CTU with the AUSM+ 
solver (top) and the RK3 with HLL solver and third-order CENO interpolation (bottom). 
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Fig. 4. — Axial Mach numbers plotted at t = 60, 120, 230 (for Muscl Hancock) and at t = 60, 120, 245 (for RK3) using dotted, 
dashed and dash-dotted lines, respectively. Because of the small-amplitude oscillations around the equilibrium position, the 
last time is not the same. The solid line gives the empirical relation l|45j l for < z < z m where z m is given by Eq. Q44I . 
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In both cases we see that the Mach number profiles plotted in Fig [4] exhibit an excellent agreement with 
the empirical relation (j4"S")) . It should be emphasized that the saturated position of the Mach disk oscillates 
around the experimental value, Eq. (|44[) . by a few percents. The oscillations are caused by disturbances 
originated at the intersection between the Mach disk and the barrel shock surrounding the jet. 



5.3. Rotated MHD Shock Tube Problem 



This test sets up a Riemann problem between left and right states given, respectively, by 
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(47) 



where parallel (||) and perpendicular (_L) components refer to the direction of structure propagation. The 
ratio of specific heats is T = 5/3. The line 2y = (A — x) (corresponding to a rotation angle of a = tan -1 2 
with respect to the y-axis) is chosen as the surface of discontinuity separating the two consta nt states and 
the computational domain [0,1] x [0,1 /iV x ] is discretized with N x x 2 zones, as in (|2005h . This setup 
yields a mesh spacing ratio Ax/ Ay = 2 thus ensuring that the initial magnetic field is divergence- free (CT 
is used to advect magnetic fields). At the boundaries in the y direction we impose translational invariance 
by applying shifted boundary condition, i.e., for any pair of indexes spanning the ghost zones, we set 
Vi y j — V^_ij+i at the lower y-boundary and Vij — V^+ij-i at the upper y-boundary. 

The structure involves three-dimensional vector fields and the outcoming wave pattern is bounded by 
two fast shocks (located x k, 0.3 and x « 0.95) enclosing two rotational discontinuities (visible only in the 
magnetic field), two slow shocks and one contact wave in at x rs 0.62. Fig. [5] shows a cut along the x 
axis at t = 0.2 cos a together with the high-resolution reference solution (solid line) computed from a one- 
dimensional integration at t — 0.2. The second order unsplit Runge-Kutta method with a CFL number of 
0.4 together with 2 nd order linear interpolation on primitive variabels and the Riemann solver of Roe were 
used for integration. 

Computations were repeated at different resolutions starting from iV x = 50 (lowest) up to N x = 1600 
(highest), by doubling the number of zones each time. The performances and accuracies of the Roe, HLLD 
and HLL Riemann solvers were compared in terms of errors (in L\ norm) and CPU time, as shown in Fig[6] 
As indicated, HLLD and Roe yields comparable errors (9.53 • 10~ 4 and 9.36 • 10" 4 at the highest resolution), 
whereas the HLL solver results in errors which are ~ 22% larger. In terms of CPU time, the HLLD solver 
is considerably faster than the Roe scheme and this is reflected in the right panel of Fig [Sj Indeed, for a 
given accuracy, the computing time offered by HLLD significantly improves over the Roe solver and it is still 
slightly better than the cheaper HLL scheme. In this sense, HLLD offers the best trade off between accuracy 
and efficiency. 

A caveat: the previous considerations hold, strictly speaking, for the particular combinations of algo- 
rithms applied to the problem in question and should not be trusted as general statements. One may find, 
for example, that if interpolation is carried on characteristic variables (instead of primitive) the overall trend 
which favours HLL over Roe (in terms of CPU time) is reversed. Still, our results illustrate how the choice 
of an algorithm or another may considerably bear upon important issues of accuracy and efficiency. 
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Fig. 5. — Rotated MHD shock-tube problem at t = 0.2 cos a. Symbols give the solution computed on a two-dimensional grid 
with 400 X 2 zone, while the solid line gives a reference solution computed for the non-rotated version at a resolution of 8192 
zones. 
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5.4. MHD Fast Rotor Problem 

A rapidly spinning (u) — 20) cylinder with higher density is embedded in a static background medium 
with uniform pressure (p = 1), threaded by a constant magnetic field (B x — 5/a/47t). Hydrodynamical 
variables smoothly change their values between the disk and the external environment using a taper function: 

{(10, lut) for r < r , 

(1 + 9/, fwr Q ) for r„ < r < n , (48) 
(1,0) otherwise. 

The disk has radius tq — 0.1 and the taper function / = {r\ — r)/{r\ — ro) terminates at r\ = 0.115. 
The initial radial velocity is zero everywhere a nd the adiabatic i ndex is T = 1.4. This problem has b een 



consid ered by several authors, see for example iBalsara k, Spicerl ( 19991 ); lLondrillo &: Del Zannal ( 20041 ); |L 



(|2005l ). and references therein. 



Because of the geometrical setting, we have setup the problem in both Cartesian and polar coordinates. 
On the Cartesian grid, the domain is the square [—5, \] 2 with outflow boundary conditions applied every- 
where. On the polar grid, we choose 0.05 < r < 0.5, < (j) < 2tt with periodic boundary conditions in cf> and 
zero-gradient at the outer radial boundary. At the inner boundary we set d r (v r /r) — d r (v^/r) = 0. For each 
geometry, we adopt the combination of algorithms listed in Table [3] and compare two different strategies to 
control the solenoidal constraint V • B = 0, i.e., the CT and Powell's methods. 

Fig. [7] shows density and magnetic pressure contours computed at t — 0.15 on the Cartesian and 
polar grids at a resolution of 400 2 zones. As the disk rotates, strong torsional Alfven waves form and 
propagate outward carrying angular momentum from the disk to the ambient. Our results fully agree 
with those previously recovered by the above-mentioned authors and computations obtained with the CT 
and Powell's scheme shows excellent agreement. Despite the higher complexity of the CT scheme where 
both staggered and zone-centered magnetic field are evolved in time, the computational cost turned out 
to be comparable (CT/Powell ~ 1.05) for each system of coordinates. Moreover, the solutions computed 
in the two different geometries show nearly identical patterns, although polar coordinates (even with their 
intrinsically higher numerical viscosity at higher radii, due to the diverging nature of the grid) better fit the 
geometrical configuration of the problem. However, computations carried out on the polar grid were ~ 5 
times slower than the Cartesian one, because of the more severe time step limitation at the inner boundary, 
where the grid narrows. These considerations can considerably affect the choice of geometry, specially in 
long-lived simulations. 



5.5. Magnetized Accretion Torus 

We now discuss an application of the code to a problem of astrophysical interest, along the lines of 



Hawlevl (l2000l) . The problem involves a magnetized, constant angular momentum (Q cx (rsin( 



in a (properly normalized) pseudo-Newtonian gravitational potential, $ = — (r — 1) . The torus has an 
equilibrium configuration described by the integral relation 



l2 



TP =C-*- 1 ,-^ T -- (49) 



(r - l)p 2 r 2 sin 2 I 

Here C is determined by the radial location of the inner edge of the torus (r = 3 in our case) and ^kep = 
r 3 / 2 /(r _ 1) ig t ne Keplerian specific angular momentum at the pressure maximum (r = 4.7) where the 
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orbital period is AT = 50. Density and pressure are tied by the polytropic relation p = Kp T , with T = 5/3. 
A poloidal magnetic field is initialized inside the torus using the <fi component of the potential vector oc 
min (p(r, 4>) — 5, 0), normalized by the condition min (2p/|B| 2 ) = 10 2 



We solve the problem in spherical as well as cylindrical coordinates. The spherical grid used for the 
problem is 1.5 < r < 20 and < 8 < ir, with 320 x 256 uniform zones covering the region 1.5 < r < 11.5 and 
7r < A8 < 37r. The remaining regions were covered by a geometrically stretched grid with 72 zones. Outflow 
boundaries are imposed at the innermost and outermost radii, whereas axisymmetry holds at 9 — and 
8 = it. In cylindrical coordinates, the computational box is < r < 20 and —20 < z < 20, with a uniform 
grid 320 x 320 in the region 1.5 < r < 11.5 and —5 < z < 5. Axisymmetry boundary conditions are imposed 
at the axis, while all the rest are set to outflow. In addition, each case is further investigated by adopting 1) 
the HLLD solver with second order slope limited reconstruction -case a and c in Table [3]- and 2) the HLL 
solver combined with piecewise parabolic reconstruction (case b and d in Table [3]). Second and third order 
Runge Kutta schemes are used for integration, respectively. 

Figure [S] shows the density logarithm at t/AT — 5.6 for the four cases respectively. Initially, the strong 
shear generates a toroidal magnetic field component leading to a pressure-driven expansion of the torus and 
accretion takes place at t/AT w 1. The growth of the pol oidal field develops into a highly turbulent regime 



accounted for by the magneto rotational instability (MRI. iBalbus fc rlawlevlll99lf ) at t/AT « 4. The effect 
is a net transport and redistribution of angular momentum. 

Even though linear reconstruction presumably introduces more diffusion than parabolic interpolation, 
the HLLD solver, being more accurate, results in computations with a higher level of fine structure (compared 
to HLL) around the tori during the turbulent phase. Besides, by normalizing the CPU time to case c), we 
found the ratio a:b and c:d to be 1 : 1.75 (in cylindrical coordinates) and 1.55 : 2.42 (in spherical coordinates), 
respectively. Indeed, since higher than 2 nd order interpolants are computationally more expensive and 
are used in conjunction with 3 rd order time accurate schemes, the HLL integration comes at an extra 
computational cost. This suggests that the use of a more accurate Riemann solver can balance the lower 
order of reconstruction in space and result, at the same time, in more efficient computations. Moreover, note 
that the equatorial symmetry breaks down more rapidly when parabolic reconstruction is employed, but it 
is retained to a higher level for the 2 nd order case. This is likely due to the larger number of operations 
introduced by the PPM reconstruction step, inevitably leading to faster-growing round off errors. 

Finally we note that computations in spherical geometry suffer from the major drawbacks of introducing 
excessive numerical viscosity at the outer radii and to severely limiting the time step at the inner boundary. 
Both aspects are indeed confirmed in our numerical tests. On the other hand, a cylindrical system of 
coordinate does not experience a similar behavior and thus a higher turbulence level is still observed at 
larger radii. 



5.6. Three-dimensional Spherically Symmetric Relativistic Shock Tube 

An initial spherical discontinuity separates a high pressure region {p = 10 3 for r < 0.4) from a uniform 
medium with p = 1 (for r > 0.4). The density is set to unity everywhere and the gas is initially static. 

Computations proceed by investigating the influence of two different equations of state (EoS), the 
constant T law with T = 5/3 (cq. [30]) and the TM equation ([31]) . Moreover, we compare the performances of 
the PPM scheme (with the Two-Shock Riemann solver) with the simpler MH (with the HLL solver) which 
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avoids characteristic decomposition of the equations (schemes a and b in Table [3| . Spherical symmetry is 
assumed in the one-dimensional calculation. 

The ensuing wave pattern (solid line in Fig [5]) is comprised of a left rarefaction wave, a contact discon- 
tinuity in the middle and a right going shock. The overall morphology is significantly affected by the choice 
of the EoS. Indeed, when the more realistic EoS is adopted, waves propagate at smaller velocities and this 
is particularly evident at the head and the tail points of the left-going rarefaction wave. This results in a 
stronger adiabatic expansion for the constant T gas and in a denser shell between the shock and the trailing 
contact wave for the TM EoS. 

Fig. [10] reports the L\ norm errors (for density) and the normalized CPU time as functions of the 
resolution Ar" 1 =50-2" (0 < n < 5) for both schemes and for the selected equations of state. Errors are 
computed against a reference solution of 16 • 10 3 zones. The cheaper scheme b) is certainly faster, although 
almost twice the resolution should be employed to achieve comparable accuracy with scheme a) . The opposite 
trend is found for the PPM scheme which is slower by almost a factor of 2 when compared to MH/HLL. 
Employment of the TM EoS results in an additional computational cost of about 25 -j- 30% independently 
of the Riemann solver. However, the simple results obtained here suggest that a more realistic EoS can 
significantly impact the solution and thus justify the direct use of Eq. (13"T|) inst ead of ([501) . Furt her reading on 



(120071) 



this to pic, together with its extension to RMHD, may be found in Ryu et al. ( 2006h : Mignone fc Mc Kinney 



Spherically symmetric problems can serve as useful verification benchmarks to tests the code perfor- 
mance on three dimensional Cartesian grids. For this purpose, we us e the second-order ChTr scheme with 
the recently developed HLLC Riemann solver ([Mignone fc Bodoll2005h (scheme c) and adopt a resolution of 



128 3 computational zones. By exploiting the symmetry, computations are carried out on one octant ([0, l] 3 ) 
only by adopting reflecting boundary conditions. Results are over-plotted using symbols in Fig. (9[ Small 
propagation structures such as the thin density shell pose serious computational challenges. Indeed we see 
that the maximum density peaks achieve ~ 85% and ~ 83% of the reference values for T = 5/3 and the TM 
equation of state, respectively. The smooth rarefaction wave is correctly captured at constant entropy with 
very small overshoots at the back. Results at higher resolution (not shown here) show increasingly better 
agreement. 

Scaling performances have been tested for this problem on different parallel platforms, all yielding 
satisfactory scaling properties. Fig. fTS] shows, for instance, the normalized execution time on the IBM SP 
Cluster 1600 p5-575 with N p — 2 P processors (1 < p < 5) for the same problem at a resolution of 128 3 , 
indicating an almost ideal scaling (oc 1/N p ). 



5.7. Three-dimensional Relativistic Magnetized Blast Wave 

A spherical region with density p = 10~ 2 and thermal pressure p = 1 is embedded in a static uniform 
medium with p = 10~ 4 and p = 3 • 10~ 5 . The sphere is centered around the origin and has radius r = 0.8. 
A linear smoothing function is applied for 0.8 < r < 1. The whole region is threaded by a constant vertical 
magnetic field in the z direction, B z = 1. The ideal equation of state with specifi c heat ratio T = 4 /3 is 



used. A similar setup in two dimensional planar geometry was earlier considered bv iKomissarovl |1999). 

Our computational grid employs 512 3 zones on the box [— 6,6] 3 . By exploiting the symmetric nature 
of the problem we reduce the computations to one octant (0 < x,y,z < 6) at half the resolution This is 
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easily achieved by symmetrizing all quantities with respect to the coordinate planes and flipping the signs 
of (B K , v K ) at y = 0, [B y , v y ) at x = 0, (v z , B x , B y ) at z = 0. 

Figs. [11] shows the magnetic field distribution, thermal pressure and Lorentz factor for the two configu- 
rations at t = 4. As a reference solution, computations on a cylindrical axisymmetric grid with [512 x 1024] 
zones are also shown. The expanding region is delimited by a fast forward shock propagating radially at 
almost the speed of light. Because of the strong magnetic confinement a jet-like structure develops in the z 
direction with a maximum Lorentz factor of 7 max ~ 4.5. This problem is particularly challenging because 
of the very small plasma (3 = 2p/\b m \ 2 = 6 • 10 -5 and instructive in checking the robustness of the code and 
the algorithm response to different kinds of degeneracies. 

Scaling tests have been conducted in a way similar to £ 15.61 and corresponding performance are shown in 
Fig. m 



6. Conclusions 

In the context of astrophysical problems, where high Mach number compressible flow dynamics plays 
a crucial role, a sizable number of numerical codes is now available to the community of scientists. Most 
of them are based on the reconstruct-solve-average strategy typical of the high-resolution shock-capturing 
Godunov-type codes. The reasons for choosing one code or another can be the most diverse, however, one 
can attempt to list a number of characteristics that a code should hold to be appealing for the skilled user, 
but not necessarily expert in numerics: 

1. a modular structure: the possibility to easily code and combine different algorithms to treat different 
physics; 

PLUTO offers a number of features which can be combined together. Besides four physical modules 
(e.g. Newtonian/relativistic hydrodynamics with or without magnetic fields), gravity, radiative cooling, 
resistivity, multidimensional geometries, equations of state may be included where required; 

2. provide the user with a number of numerical schemes tested against the most severe benchmarks; 
several algorithms (e.g. Riemann solvers, interpolations, choice of boundary conditions and so on) have 
been coded in PLUTO and their choice is dictated by the problem at hand and/or by efficiency and 
robustness issues. 

3. portability; 

we have successfully ported PLUTO to a number of different Unix-based systems among which IBM 
sp5/sp4, SGI Irix, Linux Beowulf clusters, Power Macintosh. In addition the code can run in either 
serial, single-processor or parallel machines. 

4. grid adaptivity: the ability to resolve flow features with different spatial and temporal scales on the same 
computational domain; 

the code provides a one-dimensional AMR integrator and multidimensional extensions are being de- 



veloped by taking advantage of the CHOMBO library ( [http://seesar.l bl.gov/ANAG/chombo/ 1- This 



feature will be distributed with future versions of the code. 



5. 



last but not least, user friendliness; 

a simple user-interface based on the Python scripting language is available to setup a physical problem 



-27- 



in a quick and self-explanatory way. The interface is conceived to minimize the coding efforts left to 
the user. 

The code together with its documentation is freely distributed and it is available at the web site 
http://plutocode.to.astro.i"t[ 

We would like to thank S. Ritta for his useful contributions to the under-expanded jet test and E. 
Beltritti for extensively testing the code performances on several platforms. The present work was partially 
supported by the European Community Marie Curie Actions - Human resource and mobility within the 
JETSET network under contract MRTN-CT-2004 005592. Numerical calculations were partly performed in 
CINECA (Bologna, Italy) thanks to the support by INAF. 

A. M. would like to thank R. Rosner, T. Linde and T. Plewa and all the people at the FLASH center at 
the University of Chicago for their helpful suggestions and discussions during the early development of the 
code. 
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N x CPU Time 

Fig. 6. — Errors in L\ norm (ei, 1 (q) = ^2 ~ it^l^XiAyj where qij and are the numerical and analytical solutions on 
the mesh) for the rotated MHD shock-tube problem at t = 0.2 cos a. In order to consider all discontinuities e was computed 
as the arithmetic average of £Li(p) and ez J1 (B z ). On the left, errors are plotted as function of the mesh size, whereas the left 
panel gives the errors as function of the CPU time, normalized to the fastest integration. Solid, dashed and dotted lines refer 
to computations carried with the Roc, HLLD and HLL Riemann solvers, respectively. 
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Fig. 7. — Contour maps for the MHD rotor problem at t = 0.15. Left and right panels show, respectively, computations 
carried with the CT and Powell's method in Cartesian coordinates (1 st and 3 rd rows from top) and polar coordinates (2 nd 
and 4 th rows from top). Thirty equally spaced contour levels are used for density ((0.483 < p < 13.21, 1 st and 2 nd rows) and 
magnetic pressure (0.0177 < |B| 2 /2 < 2.642, 3 rd and 4 th ros). 
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Fig. 8. — Density logarithm for the magnetized accretion torus at t/At = 5.6. Computations with the HLLD Riemann solver 
(left) and the HLL solver (right) are shown in spherical (top) and cylindrical (bottom) geometries. The region < r < 1.5 is 
excluded from computations. 
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Fig. 9. — Spherically symmetric relativistic shock tube at t = 0.4. Solid and dashed lines give the one dimensional reference 
solutions for the TM and Y = 5/3 equations of state, respectively. Density (p), pressure logarithm (logp), radial velocity (v r ) 
and entropy logarithm are shown. Over-plotted filled circles and empty boxes show the same quantities along the main diagonal 
of the three-dimensional simulation. The CFL number is 0.8 and arc 128 3 have been used. 
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Fig. 10. — Li norm errors (on the left) and CPU time (on the right) for the one dimensional spherically symmetric relativistic 
shock tube as function of the resolution. Solid and dotted lines label computations carried with the PPM scheme (a) and the 
MH scheme (b). Results pertaining to different equations of state are marked with filled circles (T = 5/3) and boxes (TM EoS). 
The CPU time has been normalized to the fastest computation. 
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Fig. 11. — Relativistic magnetized blast wave at t = 0.4 in Cartesian (left panels, cuts in the y — z plane at x = 0) and 
cylindrical (right panels) coordinates. Shown in the uppermost panels are the cylindrical radial and vertical magnetic field 
distributions with equally spaced color levels ranging from —0.052 to 0.052 (for B r ) and from 0.74 to 1.17 (for B z ). The bottom 
panels show, respectively, thermal pressure (in log scale) and the Lorentz factor with contours ranging from —4.54 to —0.74 
and from 1 to 4.45. Magnetic field lines are over-plotted on the Lorentz factor distribution. 
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Fig. 12. — Scaling of PLUTO on 2, 4, 8, 16 and 32 processors for the 3D relativistic shock tube interaction (stars) and the 
RMHD blast wave (triangles) test problems, respectively normalized to the 32- and 16- processor execution time. A grid of 128 
was used. The perfect scaling slope (<x l/N p where N p is the number of processors) is shown as the solid line. Computation 
were performed on the IBM SP Cluster 1600 p5-575. 



